﻿#include "RBF.h"

using namespace std;

constexpr double M_PI = 3.14159265358979323846;

dvector shift = dvector(0, 0, 0);
double lamda[299][300], kesai[299][300];
double* w, * x;


void printCurrentTime(string printinfo) {
    // 获取当前系统时间
    std::time_t now = std::time(nullptr);

    // 将时间转换为本地时间（线程安全）
    std::tm timeinfo;
    localtime_s(&timeinfo, &now);

    // 打印当前时间
    std::cout << printinfo << ": "
        << (timeinfo.tm_year + 1900) << "-"  // 年份
        << (timeinfo.tm_mon + 1) << "-"     // 月份
        << timeinfo.tm_mday << " "          // 日期
        << timeinfo.tm_hour << ":"          // 小时
        << timeinfo.tm_min << ":"           // 分钟
        << timeinfo.tm_sec                  // 秒
        << std::endl;
}

void generate_Gaussian_nodes()
{

    x = new double[300];

    w = new double[300];

    int i, j;

    for (i = 0; i < 298; i++)
    {

        gaulen(-1., 1., x, w, i + 2);

        for (j = 0; j <= i + 1; j++)
        {

            lamda[i][j] = w[j + 1];

            kesai[i][j] = x[j + 1];
        };
    }

    delete[] w;

    delete[] x;
}

double Cal_Fit_c(int Nums_1D, double length_lambda)
{
    if (Nums_1D == 3)
    {
        double Intercept = 0.18069, B1 = -7.02945, B2 = 123.45834, B3 = -956.11453, B4 = 4108.10754, B5 = -9248.47491, B6 = 8455.6736;
        return Intercept + B1 * length_lambda + B2 * pow(length_lambda, 2) + B3 * pow(length_lambda, 3) + B4 * pow(length_lambda, 4) + B5 * pow(length_lambda, 5) + B6 * pow(length_lambda, 6);
    }
    else if (Nums_1D == 4)
    {
        double Intercept = 29.95076, B1 = -396.10989, B2 = 2160.37682, B3 = -6220.24612, B4 = 10021.31918, B5 = -8571.32172, B6 = 3040.76241;
        return Intercept + B1 * length_lambda + B2 * pow(length_lambda, 2) + B3 * pow(length_lambda, 3) + B4 * pow(length_lambda, 4) + B5 * pow(length_lambda, 5) + B6 * pow(length_lambda, 6);
    }
    else if (Nums_1D == 5)
    {
        double Intercept = 232.59376, B1 = -1737.51008, B2 = 5322.59744, B3 = -8588.05217, B4 = 7721.64049, B5 = -3674.02415, B6 = 723.38973;
        return Intercept + B1 * length_lambda + B2 * pow(length_lambda, 2) + B3 * pow(length_lambda, 3) + B4 * pow(length_lambda, 4) + B5 * pow(length_lambda, 5) + B6 * pow(length_lambda, 6);
    }
    else if (Nums_1D == 6)
    {
        double Intercept = 9493.23661, B1 = -45371.35454, B2 = 90164.22251, B3 = -95377.92377, B4 = 56652.24064, B5 = -17917.10747, B6 = 2357.37688;
        return Intercept + B1 * length_lambda + B2 * pow(length_lambda, 2) + B3 * pow(length_lambda, 3) + B4 * pow(length_lambda, 4) + B5 * pow(length_lambda, 5) + B6 * pow(length_lambda, 6);
    }
    else if (Nums_1D == 7)
    {
        double Intercept = -4193.8859, B1 = 11921.89502, B2 = -13396.19048, B3 = 7336.06779, B4 = -1867.88249, B5 = 124.42432, B6 = 18.00615;
        return Intercept + B1 * length_lambda + B2 * pow(length_lambda, 2) + B3 * pow(length_lambda, 3) + B4 * pow(length_lambda, 4) + B5 * pow(length_lambda, 5) + B6 * pow(length_lambda, 6);
    }
    else if (Nums_1D == 8)
    {
        double Intercept = 42812.20981, B1 = -127853.52829, B2 = 158354.20542, B3 = -104182.02966, B4 = 38418.79526, B5 = -7532.38845, B6 = 613.60936;
        return Intercept + B1 * length_lambda + B2 * pow(length_lambda, 2) + B3 * pow(length_lambda, 3) + B4 * pow(length_lambda, 4) + B5 * pow(length_lambda, 5) + B6 * pow(length_lambda, 6);
    }
    else
    {
        std::cout << "Nums_1D is Error" << std::endl;
        return 0;
    }
}
void generate_NewChebyshev_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ)
{
    auto GenerateChebyshevInterpolationPoint = [=](double L, int N, int i) -> double {
        return (L / 2.0) * (cos((2.0 * (N - 1 - i) + 1) * M_PI / (2.0 * N)) + 1.0);
        };

    auto GenerateClenshawCurtisInterpolationPoint = [=](double L, int N, int i) -> double {
        return (L / 2.0) * (cos((N - 1 - i) * M_PI / N) + 1.0);
        };

    zeroX.clear();
    zeroY.clear();
    zeroZ.clear();

    double dl_x, dl_y, dl_z, dl_x_2, dl_y_2, dl_z_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[Nx];

    dl_x = length_x / (Nx - 1);
    dl_y = length_y / (Ny - 1);
    dl_z = length_z / (Nz - 1);

    dl_x_2 = dl_x / 2;
    dl_y_2 = dl_y / 2;
    dl_z_2 = dl_z / 2;

    for (int k = 0; k < Nx; k++)
    {
        double pointx = GenerateChebyshevInterpolationPoint(length_x, Nx, k);
        zeroX.push_back(pointx);
    }
    for (int k = 0; k < Ny; k++)
    {
        double pointy = GenerateChebyshevInterpolationPoint(length_y, Ny, k);
        zeroY.push_back(pointy);
    }
    for (int k = 0; k < Nz; k++)
    {
        double pointz = GenerateChebyshevInterpolationPoint(length_z, Nz, k);
        zeroZ.push_back(pointz);
    }
    zeroX[0] = 0.0;
    zeroY[0] = 0.0;
    zeroZ[0] = 0.0;
    zeroX[zeroX.size() - 1] = length_x;
    zeroY[zeroY.size() - 1] = length_y;
    zeroZ[zeroZ.size() - 1] = length_z;

    n = 0;
    for (int i = 0; i < Nx; i++)
    {
        for (int j = 0; j < Ny; j++)
        {
            for (int k = 0; k < Nz; k++)
            {
                r[n] = dvector(zeroX[i], zeroY[j], zeroZ[k]);
                n++;
            }
        }
    }
}

void generate_Chebyshev_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ)
{
    auto GenerateChebyshevInterpolationPoint = [=](double L, int N, int i) -> double {
        return (L / 2.0) * (cos((2.0 * (N - 1 - i) + 1) * M_PI / (2.0 * N)) + 1.0);
        };

    auto GenerateClenshawCurtisInterpolationPoint = [=](double L, int N, int i) -> double {
        return (L / 2.0) * (cos((N - 1 - i) * M_PI / N) + 1.0);
        };

    zeroX.clear();
    zeroY.clear();
    zeroZ.clear();

    double dl_x, dl_y, dl_z, dl_x_2, dl_y_2, dl_z_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[Nx];

    dl_x = length_x / (Nx - 1);
    dl_y = length_y / (Ny - 1);
    dl_z = length_z / (Nz - 1);

    dl_x_2 = dl_x / 2;
    dl_y_2 = dl_y / 2;
    dl_z_2 = dl_z / 2;

    for (int k = 0; k < Nx; k++)
    {
        double pointx = GenerateChebyshevInterpolationPoint(length_x, Nx, k);
        zeroX.push_back(pointx);
    }
    for (int k = 0; k < Ny; k++)
    {
        double pointy = GenerateChebyshevInterpolationPoint(length_y, Ny, k);
        zeroY.push_back(pointy);
    }
    for (int k = 0; k < Nz; k++)
    {
        double pointz = GenerateChebyshevInterpolationPoint(length_z, Nz, k);
        zeroZ.push_back(pointz);
    }

    n = 0;
    for (int i = 0; i < Nx; i++)
    {
        for (int j = 0; j < Ny; j++)
        {
            for (int k = 0; k < Nz; k++)
            {
                r[n] = dvector(zeroX[i], zeroY[j], zeroZ[k]);
                n++;
            }
        }
    }
}

void generate_tartan_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ)
{
    zeroX.clear();
    zeroY.clear();
    zeroZ.clear();

    double dl_x, dl_y, dl_z, dl_x_2, dl_y_2, dl_z_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[Nx];

    dl_x = length_x / (Nx - 1);
    dl_y = length_y / (Ny - 1);
    dl_z = length_z / (Nz - 1);

    dl_x_2 = dl_x / 2;
    dl_y_2 = dl_y / 2;
    dl_z_2 = dl_z / 2;

    for (int k = 0; k < Nx; k++)
    {
        zeroX.push_back(dl_x * k);
    }
    for (int k = 0; k < Ny; k++)
    {
        zeroY.push_back(dl_y * k);
    }
    for (int k = 0; k < Nz; k++)
    {
        zeroZ.push_back(dl_z * k);
    }

    n = 0;
    for (int i = 0; i < Nx; i++)
    {
        for (int j = 0; j < Ny; j++)
        {
            for (int k = 0; k < Nz; k++)
            {
                r[n] = dvector(zeroX[i], zeroY[j], zeroZ[k]);
                n++;
            }
        }
    }
}

void generate_tartan_interpolation_zeros(dvector* r, int N, double length, std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ)
{
    zeroX.clear();
    zeroY.clear();
    zeroZ.clear();

    double dl, dl_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[N];

    dl = length / (N - 1);

    dl_2 = dl / 2;

    for (k = 0; k < N; k++)
    {

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*N)*_pi))/(cos(_pi/(2.*N)))); //Li's chebyshev grid

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*N)*_pi));//chebyshev node

        // zeros[k]=length/2.*(1-cos(k*_pi/(N-1)));

        // zeros[k]=length/2.*(1+kesai[N-2][k]);  //gauss integral point

        zeros[k] = k * dl; // uniform

        zeroX.push_back(k * dl);
        zeroY.push_back(k * dl);
        zeroZ.push_back(k * dl);
    }

    n = 0;

    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            for (k = 0; k < N; k++)
            {
                // r[n]=dvector(i*dl,j*dl,k*dl);
                r[n] = dvector(zeros[i], zeros[j], zeros[k]);
                n++;
            }
        }
    }
}

void generate_interpolation_zeros(dvector* r, int Nx, int Ny, int Nz, double length_x, double length_y, double length_z,
    std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ)
{
    zeroX.clear();
    zeroY.clear();
    zeroZ.clear();

    double dl_x, dl_y, dl_z, dl_x_2, dl_y_2, dl_z_2;

    dl_x = length_x / (Nx - 1);
    dl_y = length_y / (Ny - 1);
    dl_z = length_z / (Nz - 1);

    dl_x_2 = dl_x / 2;
    dl_y_2 = dl_y / 2;
    dl_z_2 = dl_z / 2;

    for (int k = 0; k < Nx; k++)
    {
        zeroX.push_back(k * dl_x);
        if (k < Nx - 1)
        {
            zeroX.push_back(dl_x_2 + k * dl_x);
        }
    }
    for (int k = 0; k < Ny; k++)
    {
        zeroY.push_back(k * dl_y);
        if (k < Ny - 1)
        {
            zeroY.push_back(dl_y_2 + k * dl_y);
        }
    }
    for (int k = 0; k < Nz; k++)
    {
        zeroZ.push_back(k * dl_z);
        if (k < Nz - 1)
        {
            zeroZ.push_back(dl_z_2 + k * dl_z);
        }
    }

    int n = 0;
    for (int i = 0; i < Nx; i++)
    {
        for (int j = 0; j < Ny; j++)
        {
            for (int k = 0; k < Nz; k++)
            {
                r[n] = dvector(zeroX[i], zeroY[j], zeroZ[k]);
                n++;
            }
        }
    }
}

void generate_interpolation_zeros(dvector* r, int N, double length, std::vector<double>& zeroX, std::vector<double>& zeroY, std::vector<double>& zeroZ)
{
    zeroX.clear();
    zeroY.clear();
    zeroZ.clear();

    double dl, dl_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    zeros = new double[N];

    dl = length / (N - 1); // ÿ����ֵ�ڵ�֮��ļ��

    dl_2 = dl / 2; // �ڵ����һ��

    for (k = 0; k < N; k++)
    {

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*N)*_pi));  //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*N)*_pi))/(cos(_pi/(2.*N)))); //Li's chebyshev grid

        // zeros[k]=length/2.*(1-cos(k*_pi/(N-1)));

        // zeros[k]=length/2.*(1+kesai[N-2][k]);//gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*N)*_pi-_pi/4.)); //tangent nodes

        zeros[k] = k * dl; // uniform �����ֵ�ڵ�����꣬���� N �����ȷֲ��Ľڵ㣻�ڵ��ֵΪ k * dl���� 0 �� length ֮���� dl Ϊ���
        zeroX.push_back(k * dl);
        zeroY.push_back(k * dl);
        zeroZ.push_back(k * dl);
    }

    n = 0;

    for (i = 0; i < N; i++)
    {

        for (j = 0; j < N; j++)
        {

            for (k = 0; k < N; k++)
            {

                // r[n]=dvector(i*dl,j*dl,k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]); // ������ά�Ľڵ�����

                n++;
            }
        }
    }

    for (k = 0; k < (N - 1); k++)
    {

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*(N-1))*_pi)); //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*(N-1))*_pi))/(cos(_pi/(2.*(N-1)))));//Li's chebyshev grid

        // zeros[k]=length/2.*(1+kesai[N-3][k]); //gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*(N-1))*_pi-_pi/4.));    //tangent nodes

        zeros[k] = dl_2 + k * dl; // uniform
        // zeroX.push_back(dl_2 + k * dl);
        // zeroY.push_back(dl_2 + k * dl);
        // zeroZ.push_back(dl_2 + k * dl);
    }

    for (i = 0; i < (N - 1); i++)
    {
        for (j = 0; j < (N - 1); j++)
        {
            for (k = 0; k < (N - 1); k++)
            {
                // r[n]=dvector(dl_2+i*dl,dl_2+j*dl,dl_2+k*dl);
                r[n] = dvector(zeros[i], zeros[j], zeros[k]);
                n++;
            }
        }
    }
    // cout << "Final n: " << n << endl;
    // // ��ӡ��ά����
    // std::cout << std::fixed << std::setprecision(4); // ���������ʽ
    // std::cout << "Generated 3D coordinates:\n";
    // for (int idx = 0; idx < n; idx++) {
    //     std::cout << "Point " << idx + 1 << ": (" << std::setprecision(16) << r[idx].vect_x << ", " << r[idx].vect_y << ", " << r[idx].vect_z
    //         << ")\n";
    // }

    // ������̬������ڴ�
    delete[] zeros;
}

void generate_new_interpolation_zeros(dvector* r, int N, double length)
{

    double dl, dl_2;

    double* zeros, _pi = 3.1415927;

    int n, i, j, k;

    int section = 4;

    zeros = new double[2 * N];

    dl = length / (N - 1);

    dl_2 = dl / 2;

    for (k = 0; k < N; k++)

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*N)*_pi));  //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*N)*_pi))/(cos(_pi/(2.*N)))); //Li's chebyshev grid

        // zeros[k]=length/2.*(1-cos(k*_pi/(N-1)));

        // zeros[k]=length/2.*(1+kesai[N-2][k]);//gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*N)*_pi-_pi/4.)); //tangent nodes

        zeros[k] = k * dl; // uniform

    n = 0;

    for (i = 0; i < N; i++)

        for (j = 0; j < N; j++)

            for (k = 0; k < N; k++)
            {

                // r[n]=dvector(i*dl,j*dl,k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]);

                n++;
            }

    j = 0;

    for (k = 0; k < (N - 1); k++)

        // zeros[k]=length/2.*(1-cos((2.*k+1.)/(2.*(N-1))*_pi)); //chebyshev nodes

        // zeros[k]=length/2.*(1-(cos((2.*k+1.)/(2.*(N-1))*_pi))/(cos(_pi/(2.*(N-1)))));//Li's chebyshev grid

        // zeros[k]=length/2.*(1+kesai[N-3][k]); //gauss integral point

        // zeros[k]=length/2.*(1+tan((k+0.5)/(2.*(N-1))*_pi-_pi/4.));    //tangent nodes

    {

        zeros[j] = k * dl + dl / section; // uniform

        j++;

        zeros[j] = k * dl + dl * (section - 1) / section; // uniform

        j++;
    }

    for (i = 0; i < 2 * (N - 1); i++)

        for (j = 0; j < 2 * (N - 1); j++)

            for (k = 0; k < 2 * (N - 1); k++)
            {

                // r[n]=dvector(dl_2+i*dl,dl_2+j*dl,dl_2+k*dl);

                r[n] = dvector(zeros[i], zeros[j], zeros[k]);

                n++;
            }

    cout << "\nactual total interpolation point is" << n;
}

void generate_GG(int num_sp, double _k, Complex* GG, dvector* r, dvector rs)
{

    int n;

    double R;

    for (n = 0; n < num_sp; n++)
    {

        R = (r[n] - rs).abs_vector();

        GG[n] = exp(Complex(0, 1) * _k * R) / R;
    }
}

// 迭代方法计算阶乘
double factorial(int n)
{
    if (n == 0 || n == 1 || n < 0)
    {
        return 1;
    }
    double result = 1;
    for (int i = 2; i <= n; ++i)
    {
        result *= i;
    }
    return result;
}

void generate_the_Monomials_matrix_fortran_format(double* A, int N, double* spx, double* spy, double* spz, double c,
    double _k, double lambda, RBF_TYPE RBF_type)
{

    int i, j;

    double* a, Ri;

    VECTOR Ri_;

    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {

            Ri_ = VECTOR((spx[i] - spx[j]), (spy[i] - spy[j]), (spz[i] - spz[j]));

            Ri = Ri_.abs_vector();

            Ri = Ri / lambda; //+++++++++++++++++++++++++++++++++++++

            if (RBF_type == IMQ) {
                A[i * N + (j + 1)] = 1. / sqrt(Ri * Ri + c * c); //Inverse multiquadric
            }
            else if (RBF_type == Bessel) {
                int order = 1;
                if (Ri < 1E-20)
                {
                    A[i * N + (j + 1)] = 1 / std::pow(2, order) / factorial(order);
                }
                else
                {
                    // 计算第一类贝塞尔函数 J_nu(x)
                    double Jnx = boost::math::cyl_bessel_j(order, c * Ri);
                    A[i * N + (j + 1)] = Jnx / (c * Ri);
                }

                if (isnan(A[i * N + (j + 1)]))
                {
                    std::cout << "A[" << i * N + (j + 1) << "] == nan\n";
                    exit(0);
                }
            }
            else if (RBF_type == Cubic) {
                A[i * N + (j + 1)] = Ri * Ri * Ri; // Cubic
            }
            else if (RBF_type == TPS) {
                A[i * N + (j + 1)] = Ri * Ri * log(Ri); // Thin plate spline
            }
            else {
                std::cerr << "Error : RBF_type is error" << std::endl;
                std::exit(EXIT_FAILURE);
            }
        }
    }
}

void generate_the_Monomials_matrix_fortran_format(vector<double>& phi_mat, int N, double* spx, double* spy, double* spz, double c,
    double _k, double lambda, RBF_TYPE RBF_type)
{

    int i, j;

    double* a, Ri;

    VECTOR Ri_;

    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {

            Ri_ = VECTOR((spx[i] - spx[j]), (spy[i] - spy[j]), (spz[i] - spz[j]));

            Ri = Ri_.abs_vector();

            Ri = Ri / lambda; //+++++++++++++++++++++++++++++++++++++

            if (RBF_type == IMQ) {
                phi_mat[i * N + j] = 1. / sqrt(Ri * Ri + c * c); // Inverse multiquadric
            }
            else if (RBF_type == Bessel) {
                int order = 1;
                if (Ri < 1E-20)
                {
                    phi_mat[i * N + j] = 1 / std::pow(2, order) / factorial(order);
                }
                else
                {
                    // 计算第一类贝塞尔函数 J_nu(x)
                    double Jnx = boost::math::cyl_bessel_j(order, c * Ri);
                    phi_mat[i * N + j] = Jnx / (c * Ri);
                }

                if (isnan(phi_mat[i * N + j]))
                {
                    std::cout << "A[" << i * N + (j + 1) << "] == nan\n";
                    exit(0);
                }
            }
            else if (RBF_type == Cubic) {
                phi_mat[i * N + j] = Ri * Ri * Ri; // Cubic
            }
            else if (RBF_type == TPS) {
                phi_mat[i * N + j] = Ri * Ri * log(Ri); // Thin plate spline
            }
            else {
                std::cerr << "Error : RBF_type is error" << std::endl;
                std::exit(EXIT_FAILURE);
            }
        }
    }
}

void scattered_point_RBF_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, vector<vector<double>> Phi_mat_Pinv,
    int N, double c,
    double* spx, double* spy, double* spz, double _k, double lambda, RBF_TYPE RBF_type)
{
    
    int n, p, q;

    double* t, temp;
    vector<double> tt;
    vector<double> a(N);

    dvector Ri_;

    double Ri;

    for (n = 0; n < N; n++)
    {
        Ri_ = dvector(xx - spx[n], yy - spy[n], zz - spz[n]);

        Ri = Ri_.abs_vector();

        Ri = Ri / lambda; //++++++++++++++++++++++++++++++++++++++++++++

        if (RBF_type == IMQ) {
            a[n] = 1 / sqrt(Ri * Ri + c * c); //Inverse multiquadric
        }
        else if (RBF_type == Bessel) {
            int order = 1;
            if (Ri < 1E-20)
            {
                a[n] = 1 / std::pow(2, order) / factorial(order);
            }
            else
            {
                // 计算第一类贝塞尔函数 J_nu(x)
                double Jnx = boost::math::cyl_bessel_j(order, c * Ri);
                a[n] = Jnx / (c * Ri);
            }

            if (isnan(a[n]))
            {
                std::cout << "a[n] == nan\n";
                exit(0);
            }
        }
        else if (RBF_type == Cubic) {
            a[n] = Ri * Ri * Ri; // Cubic
        }
        else if (RBF_type == TPS) {
            a[n] = Ri * Ri * log(Ri); // Thin plate spline
        }
        else {
            std::cerr << "Error : RBF_type is error" << std::endl;
            std::exit(EXIT_FAILURE);
        }
    }

    Ga = COMPLEX(0);

    for (p = 0; p < N; p++)
    {

        tt = Phi_mat_Pinv[p];

        temp = 0;

        for (q = 0; q < N; q++)
        {

            temp = temp + tt[q] * a[q];
        }

        Ga = Ga + temp * GG[p];
    }
}

void scattered_point_RBF_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, double* P[],
    double* a, int N, double c,
    double* spx, double* spy, double* spz, double _k, double lambda, RBF_TYPE RBF_type)
{

    int n, p, q;

    double* t, temp;

    dvector Ri_;

    double Ri;

    for (n = 0; n < N; n++)
    {

        Ri_ = dvector(xx - spx[n], yy - spy[n], zz - spz[n]);

        Ri = Ri_.abs_vector();

        Ri = Ri / lambda; //++++++++++++++++++++++++++++++++++++++++++++

        // if(Ri<1e-20)

        //  a[n]=0;

        // else

        //  a[n]=Ri*Ri*log(Ri);  //Thin plate spline

        // if(Ri<1e-20)

        //   a[n]=0;

        // else

        //	  a[n]=Ri*Ri*Ri;  //Cubic

        // a[n]=pow(sqrt(c*c+Ri*Ri),3);  //Natural cubic spline

        // a[n]=log(c*c+Ri*Ri); //Multilog

        //  a[n]=Ri;  //Linear

        //  a[n]=1/sqrt(Ri*Ri+c*c); //Inverse multiquadric

        // a[n]=cos(_k*Ri)/sqrt(Ri*Ri+c*c);

        // a[n]=sqrt(Ri*Ri+c*c); //Multiquadric

        // a[n]=exp(-c*Ri*Ri); //Gaussian

        // a[n]=1/(Ri*Ri+c*c);//Inverse quadratic

        // double alpha = 0.4;
        // a[n] = (1 - alpha) * exp(-c*c*Ri*Ri) + alpha * (1.0/sqrt(Ri*Ri+c*c));

        if (RBF_type == IMQ) {
            a[n] = 1 / sqrt(Ri * Ri + c * c); //Inverse multiquadric
        }
        else if (RBF_type == Bessel) {
            int order = 1;
            if (Ri < 1E-20)
            {
                a[n] = 1 / std::pow(2, order) / factorial(order);
            }
            else
            {
                // 计算第一类贝塞尔函数 J_nu(x)
                double Jnx = boost::math::cyl_bessel_j(order, c * Ri);
                a[n] = Jnx / (c * Ri);
            }

            if (isnan(a[n]))
            {
                std::cout << "a[n] == nan\n";
                exit(0);
            }
        }
        else if (RBF_type == Cubic) {
            a[n] = Ri * Ri * Ri; // Cubic
        }
        else if (RBF_type == TPS) {
            a[n] = Ri * Ri * log(Ri); // Thin plate spline
        }
        else {
            std::cerr << "Error : RBF_type is error" << std::endl;
            std::exit(EXIT_FAILURE);
        }

        // if (Ri > 1) // with compact support

        //    a[n] = 0;

        // else

        //    // a[n]=pow((1-Ri),4)*(4*Ri+1); //C2

        //    a[n] = pow((1 - Ri), 6) * (35 * Ri * Ri + 18 * Ri + 3); // C4
    }

    Ga = COMPLEX(0);

    for (p = 0; p < N; p++)
    {

        t = P[p];

        temp = 0;

        for (q = 0; q < N; q++)
        {

            temp = temp + t[q] * a[q];
        }

        Ga = Ga + temp * GG[p];
    }
}

double LagrangeBaseFunction(int idx, int k, double* tk, double t) // lagrange base function
{
    double ret = 1.0;

    for (int i = 0; i < k; i++)
    {
        if (idx != i)
            ret *= (t - tk[i]) / (tk[idx] - tk[i]);
    }

    return ret;
}

void scattered_point_lagrange_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG,
    int Nums_1D_x, int Nums_1D_y, int Nums_1D_z,
    vector<double> zeroX, vector<double> zeroY, vector<double> zeroZ)
{
    for (int i = 0; i < Nums_1D_x; i++)
    {
        for (int j = 0; j < Nums_1D_y; j++)
        {
            for (int k = 0; k < Nums_1D_z; k++)
            {
                int q = i * Nums_1D_y * Nums_1D_z + j * Nums_1D_z + k;
                double temp = LagrangeBaseFunction(i, Nums_1D_x, zeroX.data(), xx) * LagrangeBaseFunction(j, Nums_1D_y, zeroY.data(), yy) * LagrangeBaseFunction(k, Nums_1D_z, zeroZ.data(), zz);
                Ga += temp * GG[q];
            }
        }
    }
}

void scattered_point_lagrange_interpolation(Complex& Ga, double xx, double yy, double zz, Complex* GG, int numsp_1D,
    vector<double> zeroX, vector<double> zeroY, vector<double> zeroZ)
{
    for (int i = 0; i < numsp_1D; i++)
    {
        for (int j = 0; j < numsp_1D; j++)
        {
            for (int k = 0; k < numsp_1D; k++)
            {
                int q = i * numsp_1D * numsp_1D + j * numsp_1D + k;
                double temp = LagrangeBaseFunction(i, numsp_1D, zeroX.data(), xx) * LagrangeBaseFunction(j, numsp_1D, zeroY.data(), yy) * LagrangeBaseFunction(k, numsp_1D, zeroZ.data(), zz);
                Ga += temp * GG[q];
            }
        }
    }
}

//使用Eigen计算伪逆，phi_mat按照行主序储存，return按照行主序储存的伪逆
vector<vector<double>> Eigen_Cal_PseudoInverse(vector<double> phi_mat, int K)
{
    // 使用 Map 将 std::vector<double> 转换为 Eigen 矩阵
    Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> Phi_Mat(phi_mat.data(), K, K);

    // 使用 BDCSVD 计算伪逆
    Eigen::MatrixXd pinv = Phi_Mat.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Eigen::MatrixXd::Identity(K, K));

    // // 使用 jacobiSVD 计算伪逆
    // Eigen::MatrixXd pinv = Phi_Mat.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Eigen::MatrixXd::Identity(K, K));

    // 按照行主序将伪逆矩阵的数据存储到 vector<vector<double>> 中
    std::vector<std::vector<double>> Phi_Mat_Pinv(K, std::vector<double>(K));
    for (int i = 0; i < K; ++i) {
        for (int j = 0; j < K; ++j) {
            Phi_Mat_Pinv[i][j] = pinv(i, j);
        }
    }
    return Phi_Mat_Pinv;
}

pair<double, double> New_test_scattered_RBF(const double length_x, const double length_y, const double length_z, const double lambda, int grid_type,
    int Nums_1D_x, int Nums_1D_y, int Nums_1D_z, const double c_new,
    vector<std::pair<vector<int>, vector<double>>>& Nums1Dxyz_Lengthxyz_OptimalC_list, RBF_TYPE RBF_type)
{
    double pi = 3.1415926535897932384626433832795028841971693993751;

    int numsp_1D; // numsp_1D: the number of scattered points in 1 dimension.

    const double miu0 = 1.2566370614359172953850573533118e-6;
    const double epsilon0 = 8.85418781762039e-12;
    const double c_wave_speed = 1.0 / sqrt(miu0 * epsilon0);

    std::cout << "X_length = " << length_x << "  Y_length = " << length_y << "  Z_length = " << length_z << "\n";
    std::cout << "X_length(/lambda) = " << length_x / lambda << "  Y_length(/lambda) = " << length_y / lambda << "  Z_length(/lambda) = " << length_z / lambda << "\n";
    std::cout << "lambda = " << lambda << std::endl;

    double c; // c: The shape parameter used in radial basis functions.

    int i, j, k;

    ///////////////////////////////////////

    // Generate the scattered points;

    ///////////////////////////////////////

    int test_type;

    test_type = grid_type;

    double RBF_min_error_qthPoint = std::numeric_limits<double>::max();
    double LBF_error_qthPoint = std::numeric_limits<double>::max();
    double RBF_min_Average_error = std::numeric_limits<double>::max();
    double RBF_optimal_c = 0.0;
    double RBF_optimal_c_qth = 0.0;
    Complex RBF_optimal_G = Complex(0.0, 0.0);
    bool Print_Interpolation_Point = false;
    double SVD_min_error = std::numeric_limits<double>::max();
    double SVD_max_error = std::numeric_limits<double>::min();
    vector<double> RBF_Min_error_list, LBF_error_list;
    double LBF_Average_error = -1;
    bool CalLBF = true;

    ReUesData reusedata;

    int num_sp; // the total number of interpolation points of a tartan grid or a staggered one.

    if (test_type == TartanGrids || test_type == ChebyshevGrids || test_type == NewChebyshevGrids)
        num_sp = Nums_1D_x * Nums_1D_y * Nums_1D_z;
    else if (test_type == StaggeredTartanGrids)
        num_sp = Nums_1D_x * Nums_1D_y * Nums_1D_z + (Nums_1D_x - 1) * (Nums_1D_y - 1) * (Nums_1D_z - 1);
    else {
        std::cout << "test_type is Errro" << std::endl;
    }

    cout << "\nnum_sp=" << num_sp << std::endl;

    // generate_Gaussian_nodes();

    dvector* rr; // Array used to store the interpolation points.

    double* spx, * spy, * spz;

    double* P[10000];

    double dc = 0.1;

    int nc = 0;

    Complex* GG, * GGx, * GGy, * GGz;

    GG = new Complex[num_sp]; //*10
    GGx = new Complex[num_sp];
    GGy = new Complex[num_sp];
    GGz = new Complex[num_sp];

    double* SA;

    double* svdA, * svdV, * svdW; //,*P; //svdA=svdU*svdW*svdVT, after svd decomposition svdA is replaced by svdU

    SA = new double[num_sp * num_sp + 1];
    for (int i = 0; i < num_sp * num_sp + 1; i++)
        SA[i] = 0.0;

    svdA = new double[(num_sp + 1) * (num_sp + 1)];
    for (int i = 0; i < (num_sp + 1) * (num_sp + 1); i++)
        svdA[i] = 0.0;

    svdV = new double[(num_sp + 1) * (num_sp + 1)];
    for (int i = 0; i < (num_sp + 1) * (num_sp + 1); i++)
        svdV[i] = 0.0;

    svdW = new double[(num_sp + 1)];
    for (int i = 0; i < (num_sp + 1); i++)
        svdW[i] = 0.0;

    double* Recovered_A = new double[num_sp * num_sp];
    double* A_new = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
    {
        Recovered_A[i] = 0.0;
        A_new[i] = 0.0;
    }
    double* svd_U = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
        svd_U[i] = 0.0;
    double* svd_S = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
        svd_S[i] = 0.0;
    double* svd_V = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
        svd_V[i] = 0.0;

    for (i = 0; i < num_sp; i++)
        P[i] = new double[num_sp]; // Allocate memory for matrix A.

    rr = new dvector[num_sp]; // *10

    spx = new double[num_sp];

    spy = new double[num_sp];

    spz = new double[num_sp];

    std::vector<double> zeroX, zeroY, zeroZ;

    std::vector<vector<double>> TestXYZ_list;

    double* poly;

    poly = new double[num_sp];

    double dlength;

    dlength = 0.1;

    int NC;

    if (c_new < 0.0)
    {
        NC = 3000;
    }
    else
    {
        NC = 1;
    }

    int Error_increase_count = 0;
    double Previous_RBF_Average_Error = -1;
    for (int nc = 0; nc < NC; nc++)
    {
        c = 0.01 + 0.01 * nc;
        if (c_new > 0.0)
        {
            c = c_new;
        }

        if (nc == 0) printCurrentTime("***$$$*** Build interpolation point start");
        if (test_type == StaggeredTartanGrids)
        {
            generate_interpolation_zeros(rr, Nums_1D_x, Nums_1D_y, Nums_1D_z, length_x, length_y, length_z, zeroX, zeroY, zeroZ);
        }
        else if (test_type == TartanGrids)
        {
            generate_tartan_interpolation_zeros(rr, Nums_1D_x, Nums_1D_y, Nums_1D_z, length_x, length_y, length_z, zeroX, zeroY, zeroZ);
        }
        else if (test_type == ChebyshevGrids) {
            generate_Chebyshev_interpolation_zeros(rr, Nums_1D_x, Nums_1D_y, Nums_1D_z, length_x, length_y, length_z, zeroX, zeroY, zeroZ);
        }
        else if (test_type == NewChebyshevGrids) {
            generate_NewChebyshev_interpolation_zeros(rr, Nums_1D_x, Nums_1D_y, Nums_1D_z, length_x, length_y, length_z, zeroX, zeroY, zeroZ);
        }
        if (nc == 0) printCurrentTime("***$$$*** Build interpolation point end");

        // for(int i = 0; i < zeroX.size(); i++) {
        //     std::cout << zeroX[i] << std::endl;
        // }

        for (i = 0; i < num_sp; i++)
        {
            spx[i] = rr[i].vect_x;
            spy[i] = rr[i].vect_y;
            spz[i] = rr[i].vect_z;

            reusedata.spx.push_back(rr[i].vect_x);
            reusedata.spy.push_back(rr[i].vect_y);
            reusedata.spz.push_back(rr[i].vect_z);
        }

        //////////////////////////////////////////////////////////////

        // Generate the Green's functions and store them in array GG.

        //////////////////////////////////////////////////////////////

        double _k = 2 * pi / lambda;

        // double distance = 2 * min(length_x, min(length_y, length_z));
        double distance_x = 2 * length_x;
        double distance_y = 2 * length_y;
        double distance_z = 2 * length_z;

        // distance:The distance from the original point to the nearest point of the nearest far group along x axis.

        dvector rs, rs_x, rs_y, rs_z; // Source point.

        if (test_type == StaggeredTartanGrids || test_type == TartanGrids || test_type == ChebyshevGrids || test_type == NewChebyshevGrids)
        {
            // rs = shift + dvector(distance, 0.5 * length_y, 0.5 * length_z);
            rs_x = shift + dvector(distance_x, 0.5 * length_y, 0.5 * length_z);
            rs_y = shift + dvector(0.5 * length_x, distance_y, 0.5 * length_z);
            rs_z = shift + dvector(0.5 * length_x, 0.5 * length_y, distance_z);
        }
        else
        {
            std::cout << "grid type error" << std::endl;
            exit(0);
            // rs = shift + dvector(distance, (0.5) * length_y, (0.5 + (dlength * nc)) * length_z);
        }

        if (abs(length_x - min(length_x, min(length_y, length_z))) < 1E-10) {
            rs = rs_x;
        }
        else if (abs(length_y - min(length_x, min(length_y, length_z))) < 1E-10) {
            rs = rs_y;
        }
        else if (abs(length_z - min(length_x, min(length_y, length_z))) < 1E-10) {
            rs = rs_z;
        }
        else {
            std::cerr << "Errro : rs != rs_x && rs != rs_y && rs != rs_z" << std::endl;
            exit(EXIT_FAILURE);
        }
        // std::cout << "rs = (" << rs.vect_x << ", " << rs.vect_y << ", " << rs.vect_z << ")\n"; 

        if (nc == 0) printCurrentTime("***$$$*** Build GG start");
        generate_GG(num_sp, _k, GG, rr, rs);
        if (nc == 0) printCurrentTime("***$$$*** Build GG end");
        // generate_GG(num_sp, _k, GGx, rr, rs_x);
        // generate_GG(num_sp, _k, GGy, rr, rs_y);
        // generate_GG(num_sp, _k, GGz, rr, rs_z);

        vector<double> phi_mat(num_sp * num_sp);
        if (nc == 0) printCurrentTime("***$$$*** Build phi mat start");
        // generate_the_Monomials_matrix_fortran_format(SA, num_sp, spx, spy, spz, c, _k, lambda, RBF_type);
        generate_the_Monomials_matrix_fortran_format(phi_mat, num_sp, spx, spy, spz, c, _k, lambda, RBF_type);
        if (nc == 0) printCurrentTime("***$$$*** Build phi mat end");

        if (nc == 0) printCurrentTime("***$$$*** Eigen Cal PseudoInverse start");
        vector<vector<double>> phi_mat_Pinv = Eigen_Cal_PseudoInverse(phi_mat, num_sp);
        if (nc == 0) printCurrentTime("***$$$*** Eigen Cal PseudoInverse end");

        // for (i = 0; i < num_sp; i++)
        // {
        //     for (j = 0; j < num_sp; j++)
        //     {
        //         svdA[(i + 1) * (num_sp + 1) + j + 1] = SA[j * num_sp + (i + 1)];
        //     }
        // }

        // if(nc == 0) printCurrentTime("***$$$*** SVD start");
        // svdcmp(svdA, num_sp, num_sp, svdW, svdV);
        // if(nc == 0) printCurrentTime("***$$$*** SVD end");

        // for (i = 0; i < num_sp; i++)
        // {
        //     for (j = 0; j < num_sp; j++)
        //     {
        //         P[i][j] = 0;

        //         for (k = 0; k < num_sp; k++)
        //         {
        //             P[i][j] += svdV[(i + 1) * (num_sp + 1) + (k + 1)] * svdA[(j + 1) * (num_sp + 1) + (k + 1)] /
        //                        svdW[k + 1];
        //         }
        //     }
        // }

        ///////////////////////////////////////////////////////////////////////////////

        // Test the accuracy of the interpolation in some uniformly spaced points.

        ///////////////////////////////////////////////////////////////////////////////

        int num_p_x = Nums_1D_x + 1; // num_p: the number of test points in one dimension.
        int num_p_y = Nums_1D_y + 1;
        int num_p_z = Nums_1D_z + 1;

        double dl_x = length_x / (num_p_x - 1);
        double dl_y = length_y / (num_p_y - 1);
        double dl_z = length_z / (num_p_z - 1);

        double xx, yy, zz; // Used to store the coordinates of the current points.

        Complex G0, Ga; // G0: the true value; Ga: the approximated value.

        dvector rf;

        double R, error, RBF_errorsum, LBF_errorsum;

        double errormin = 1e10, errormax = -1e10;
        vector<double> min_error_XYZ, max_error_XYZ;
        vector<double> RBF_error_list;

        RBF_errorsum = 0.0;
        LBF_errorsum = 0.0;

        TestXYZ_list.clear();
        for (i = 0; i < num_p_x; i++) {
            for (j = 0; j < num_p_y; j++) {
                for (k = 0; k < num_p_z; k++) {
                    TestXYZ_list.push_back(vector<double>{i* dl_x, j* dl_y, k* dl_z});
                }
            }
        }

        if (nc == 0) printCurrentTime("***$$$*** Test interpolation point average error start");
        if (nc == 0) std::cout << "TestXYZ_list.size() = " << TestXYZ_list.size() << std::endl;

        //Test LBF Error
        for (int i = 0; i < TestXYZ_list.size(); i++) 
        {       
            if (LBF_error_list.size() < TestXYZ_list.size())
            {
                double Test_x = TestXYZ_list[i][0];
                double Test_y = TestXYZ_list[i][1];
                double Test_z = TestXYZ_list[i][2];

                dvector rf_new = dvector(Test_x, Test_y, Test_z);

                double R = (rs - rf_new).abs_vector();
                Complex Real_G = exp(Complex(0, 1) * _k * R) / R;
                Complex LBF_G = Complex(0.0, 0.0);

                scattered_point_lagrange_interpolation(LBF_G, Test_x, Test_y, Test_z, GG, Nums_1D_x, Nums_1D_y, Nums_1D_z, zeroX, zeroY, zeroZ);
                double LBF_error = abs(Real_G - LBF_G) / abs(Real_G);

                LBF_error_list.push_back(LBF_error);
                LBF_errorsum += LBF_error;
                if (i == TestXYZ_list.size() - 1)
                {
                    LBF_Average_error = LBF_errorsum / TestXYZ_list.size();
                }
            }
            else break;
        }

#pragma omp parallel for reduction(+:RBF_errorsum) schedule(static, 10)
        for (int i = 0; i < TestXYZ_list.size(); i++)
        {
            double Test_x = TestXYZ_list[i][0];
            double Test_y = TestXYZ_list[i][1];
            double Test_z = TestXYZ_list[i][2];

            dvector rf_new = dvector(Test_x, Test_y, Test_z);

            double R = (rs - rf_new).abs_vector();
            Complex Real_G = exp(Complex(0, 1) * _k * R) / R;
            Complex RBF_G = Complex(0.0, 0.0);
            // scattered_point_RBF_interpolation(RBF_G, Test_x, Test_y, Test_z, GG, P, poly, num_sp, c, spx, spy, spz, _k, lambda, RBF_type);
            scattered_point_RBF_interpolation(RBF_G, Test_x, Test_y, Test_z, GG, phi_mat_Pinv, num_sp, c, spx, spy, spz, _k, lambda, RBF_type);
            double RBF_error = abs(Real_G - RBF_G) / abs(Real_G);

            RBF_errorsum += RBF_error;
        }

        double RBF_Average_error = RBF_errorsum / TestXYZ_list.size();


        if (RBF_min_Average_error - RBF_Average_error > 1E-10) {
            RBF_min_Average_error = RBF_Average_error;
            RBF_optimal_c = c;
        }

        if ((0.02 - RBF_Average_error) > 1E-10)
        {
            std::cout << "RBF_min_Average_error < 0.02" << std::endl;
            break;
        }

        if (Previous_RBF_Average_Error < 0.0)
            Previous_RBF_Average_Error = RBF_Average_error;
        else if (RBF_Average_error - Previous_RBF_Average_Error > 1E-10)
        {
            Error_increase_count++;
            if (Error_increase_count == 30)
            {
                std::cout << "Error_increase_count == 30" << std::endl;
                break;
            }
        }
        else
        {
            Error_increase_count = 0;
        }
        Previous_RBF_Average_Error = RBF_Average_error;
    }
    if (nc == 0) printCurrentTime("***$$$*** Test interpolation point average error end");

    if (0.02 - RBF_min_Average_error > 1E-10) {
        Nums1Dxyz_Lengthxyz_OptimalC_list.push_back(make_pair(vector<int>{Nums_1D_x, Nums_1D_y, Nums_1D_z}, vector<double>{length_x, length_y, length_z, RBF_optimal_c, RBF_min_Average_error, LBF_Average_error}));
    }

    // std::cout << "RBF_min_error_qthPoint = " << RBF_min_error_qthPoint << "   LBF_error_qthPoint = " << LBF_error_qthPoint << "\n";
    // std::cout << "RBF_optimal_c_qth = " << RBF_optimal_c_qth << "\n\n";

    // if(RBF_min_Average_error - 0.02 > 1E-10) {
    std::cout << "RBF_min_Average_error = " << std::setprecision(12) << RBF_min_Average_error << "   ";
    std::cout << "LBF_Average_error = " << std::setprecision(12) << LBF_Average_error << "\n";
    std::cout << "RBF_optimal_c = " << std::setprecision(6) << RBF_optimal_c << "\n\n";
    // for(int i = 0; i < RBF_Min_error_list.size(); i++) {
    //     std::cout << "RBF_Min_error_list[" << i << "] = " << RBF_Min_error_list[i] << "   ";
    //     std::cout << "LBF_error_list[" << i << "] = " << LBF_error_list[i] << "\n";
    // }
    // std::cout << "RBF_optimal_G = " << std::setprecision(18) << RBF_optimal_G << "\n";
    // std::cout << "\nsvd_min_error = " << std::setprecision(18) << SVD_min_error << "   svd_max_error = " << std::setprecision(18) << SVD_max_error << "\n";
    std::cout << "=============================================================\n\n";
    // }

    delete[] GG, GGx, GGy, GGz;
    delete[] SA, svdA, svd_V, svdW, Recovered_A, A_new, svd_U, svd_S, svd_V;
    delete[] rr, spx, spy, spz, poly;
    for (int i = 0; i < num_sp; i++) {
        delete[] P[i];
    }

    pair<double, double> RBF_LBF_Error = make_pair(RBF_min_Average_error, LBF_Average_error);

    return RBF_LBF_Error;
    // return LBF_Average_error;
}

double new_test_scattered_lagrange(const double length_new, const double lambda, int grid_type, int Nums_1D, const double c_new, std::vector<std::pair<double, double>>& Length_OptimalC_list)
{

    double pi = 3.1415926535897932384626433832795028841971693993751;

    int numsp_1D; // numsp_1D: the number of scattered points in 1 dimension.

    // double length=3,lambda=1;

    const double miu0 = 1.2566370614359172953850573533118e-6;
    const double epsilon0 = 8.85418781762039e-12;
    const double c_wave_speed = 1.0 / sqrt(miu0 * epsilon0);
    // double lambda = c_wave_speed / freq;
    // double lambda = 0.333333333;

    double length = length_new;

    // double lambda = 0.333333333;

    // ofstream std::cout;
    std::cout << "X_length = " << length << "  Y_length = " << length << "  Z_length = " << length << "\n";
    std::cout << "X_length(/lambda) = " << length / lambda << "  Y_length(/lambda) = " << length / lambda << "  Z_length(/lambda) = " << length / lambda << "\n";
    std::cout << "lambda = " << lambda << std::endl;

    double c; // c: The shape parameter used in radial basis functions.

    int i, j, k;

    ///////////////////////////////////////

    // Generate the scattered points;

    ///////////////////////////////////////

    int test_type;

    // test_type=1: Test the average interpolation error of the green's function generated

    // from the source point that is one cube length away from the field cube using the staggered tartan grid.

    // test_type=2: Test the average interpolation error of the green's function generated

    // from the source point that is one cube length away from the field cube using the ordinary tartan grid.

    // test_type=3: Test the average interpolation error of the green's function generated

    // from the source one cube length away from the field cube using the staggered tartan grid,

    // while the source point is moving from S to S' as in Fig. 7 of the paper.

    // test_type=4: Test the average interpolation error of the green's function generated

    // from the source one cube length away from the field cube using the ordinary tartan grid,

    // while the source point is moving from S to S' as in Fig. 7 of the paper.

    bool constant_c = true;

    // constant_c=false: the shape parameter c is varied.

    // constant_c=true: the shape parameter c is constant and should be input from the computer.

    // cout << "\ntest_type=1: Test the average interpolation error of the green's function generated"

    //     << "\nfrom the source point that is one cube length away from the field cube using the staggered tartan grid."

    //     << "\ntest_type=2: Test the average interpolation error of the green's function generated"

    //     << "\nfrom the source point that is one cube length away from the field cube using the ordinary tartan grid."

    //     << "\ntest_type=3: Test the average interpolation error of the green's function generated"

    //     << "\nfrom the source one cube length away from the field cube using the staggered tartan grid,"

    //     << "\nwhile the source point is moving from S to S' as in Fig. 7 of the paper."

    //     << "\ntest_type=4: Test the average interpolation error of the green's function generated"

    //     << "\nfrom the source one cube length away from the field cube using the ordinary tartan grid,"

    //     << "\nwhile the source point is moving from S to S' as in Fig. 7 of the paper."

    //     << "\ntest_type=5: Test the average interpolation error of the green's function generated"

    //     << "\nfrom the source point that is one cube length away from the field cube using the new tartan grid.";

    // cout << "\n test_type=";

    // cin >> test_type;

    // cout << "\nplease input 1D interpolation number";

    // cin >> numsp_1D;

    //// cout<<"\nIs the parameter c constant(Yes:1|No:0)";

    //// cin>>constant_c;

    // cout << "\nplease input parameter c";

    // cin >> c;

    test_type = grid_type;
    numsp_1D = Nums_1D;
    c = 1.0;

    double RBF_min_error_qthPoint = std::numeric_limits<double>::max();
    double LBF_error_qthPoint = std::numeric_limits<double>::max();
    double RBF_min_Average_error = std::numeric_limits<double>::max();
    double RBF_optimal_c = 0.0;
    double RBF_optimal_c_qth = 0.0;
    Complex RBF_optimal_G = Complex(0.0, 0.0);
    bool Print_Interpolation_Point = false;
    double SVD_min_error = std::numeric_limits<double>::max();
    double SVD_max_error = std::numeric_limits<double>::min();
    vector<double> RBF_Min_error_list, LBF_error_list;
    double LBF_Average_error = -1;
    bool CalLBF = true;

    c = c * length / lambda; //@@@@@@@@@@@@@@@

    // cout<<"\nThe number of intepolation point along one dimension is";

    // cin>>numsp_1D;

    int num_sp; // the total number of interpolation points of a tartan grid or a staggered one.

    if (test_type == 2 || test_type == 4)
        num_sp = numsp_1D * numsp_1D * numsp_1D;

    else if (test_type == 1 || test_type == 3)
        num_sp = numsp_1D * numsp_1D * numsp_1D + (numsp_1D - 1) * (numsp_1D - 1) * (numsp_1D - 1);

    else if (test_type == 5)
        num_sp = 8 * (numsp_1D - 1) * (numsp_1D - 1) * (numsp_1D - 1) + numsp_1D * numsp_1D * numsp_1D;

    cout << "\nnum_sp=" << num_sp << std::endl;

    generate_Gaussian_nodes();

    // cout<<"\nThe wave length is(Note the cube length is unit):";

    // cin>>lambda;

    // length: the cube length.

    // lambda: the wavelength of the EM wave.

    dvector* rr; // Array used to store the interpolation points.

    double* spx, * spy, * spz;

    double* P[10000];

    double dc = 0.1;

    int nc = 0;

    ofstream fo;

    ofstream fot; // test file

    // fot.open("testfile.txt");

    Complex* GG;

    GG = new Complex[num_sp]; //*10

    double* SA;

    double* svdA, * svdV, * svdW; //,*P; //svdA=svdU*svdW*svdVT, after svd decomposition svdA is replaced by svdU

    SA = new double[num_sp * num_sp + 1];
    for (int i = 0; i < num_sp * num_sp + 1; i++)
        SA[i] = 0.0;

    svdA = new double[(num_sp + 1) * (num_sp + 1)];
    for (int i = 0; i < (num_sp + 1) * (num_sp + 1); i++)
        svdA[i] = 0.0;

    svdV = new double[(num_sp + 1) * (num_sp + 1)];
    for (int i = 0; i < (num_sp + 1) * (num_sp + 1); i++)
        svdV[i] = 0.0;

    svdW = new double[(num_sp + 1)];
    for (int i = 0; i < (num_sp + 1); i++)
        svdW[i] = 0.0;

    double* Recovered_A = new double[num_sp * num_sp];
    double* A_new = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
    {
        Recovered_A[i] = 0.0;
        A_new[i] = 0.0;
    }
    double* svd_U = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
        svd_U[i] = 0.0;
    double* svd_S = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
        svd_S[i] = 0.0;
    double* svd_V = new double[num_sp * num_sp];
    for (int i = 0; i < num_sp * num_sp; i++)
        svd_V[i] = 0.0;

    for (i = 0; i < num_sp; i++)

        P[i] = new double[num_sp]; // Allocate memory for matrix A.

    rr = new dvector[num_sp]; // *10

    spx = new double[num_sp];

    spy = new double[num_sp];

    spz = new double[num_sp];

    std::vector<double> zeroX, zeroY, zeroZ;

    std::vector<vector<double>> TestXYZ_list;

    double* poly;

    poly = new double[num_sp];

    double dlength;

    dlength = 0.1;

    int NC;

    if (test_type == 1 || test_type == 2 || test_type == 5)
        NC = 1000;

    else
        NC = 15;

    if (c_new < 0.0)
    {
        NC = 1000;
    }
    else
    {
        NC = 1;
    }

    int Error_increase_count = 0;
    double Previous_RBF_Average_Error = -1;
    for (int nc = 599; nc < NC; nc++)
    {
        // for(; c <= max_c; c += c_step) {
        // c = c / lambda;
        // cout<<"\nnc="<<nc;

        // nc=30;

        // if(test_type==3||test_type==4||constant_c)

        //{

        //	cout<<"\n Please input a constant shape parameter C:";

        //   cin>>c;

        //}

        // else

        // c=(0.5+nc*dc)*length/(numsp_1D-1);  //c: used in radial basis functions.

        // c=length/(numsp_1D-1);

        c = 0.01 + 0.01 * nc;
        if (c_new > 0.0)
        {
            c = c_new;
        }

        if (test_type == 1 || test_type == 3)

            generate_interpolation_zeros(rr, numsp_1D, length, zeroX, zeroY, zeroZ);

        else if (test_type == 2 || test_type == 4)

            generate_tartan_interpolation_zeros(rr, numsp_1D, length, zeroX, zeroY, zeroZ);

        else if (test_type == 5)

            generate_new_interpolation_zeros(rr, numsp_1D, length);

        for (i = 0; i < num_sp; i++)
        {

            spx[i] = rr[i].vect_x;

            spy[i] = rr[i].vect_y;

            spz[i] = rr[i].vect_z;
        }

        /*if (!Print_Interpolation_Point) {
            for (i = 0; i < num_sp; i++) {
                std::cout << "x=" << spx[i] << "  y=" << spy[i] << "  z=" << spz[i] << "\n";
            }
            Print_Interpolation_Point = true;
        }*/

        //////////////////////////////////////////////////////////////

        // Generate the Green's functions and store them in array GG.

        //////////////////////////////////////////////////////////////

        double _k = 2 * pi / lambda;

        double distance = 2 * length;

        // distance:The distance from the original point to the nearest point of the nearest far group along x axis.

        dvector rs; // Source point.

        if (test_type == 1 || test_type == 2 || test_type == 5)

            rs = shift + dvector(distance, 0.5 * length, 0.5 * length);

        else
            rs = shift + dvector(distance, (0.5) * length, (0.5 + (dlength * nc)) * length);

        // std::cout << "\nrs = (" << rs.vect_x << ", " << rs.vect_y << ", " << rs.vect_z << ")\n";

        generate_GG(num_sp, _k, GG, rr, rs); // rr�洢�����еĲ�ֵ�㣬�������������rs�����еĲ�ֵ��֮������и��ֺ���Ȼ�󱣴�����GG��

        // generate_GG(_k,GG,rr,rs,numsp_1D);

        // for(i=0;i<num_sp;i++)

        //{

        //   fot<<i<<" Green="<<real(GG[i])<<endl;

        //}

        ////////////////////////////////////////

        // Generate the Monomials matrix.

        ////////////////////////////////////////

        // for(i=0;i<num_sp;i++)

        //{

        //	A[i]=new double[num_sp];  //Allocate memory for matrix A.

        //}

        // generate_the_Monomials_matrix(A,num_sp,spx,spy,spz);

        // generate_the_Monomials_matrix(A,num_sp,spx,spy,spz,c);

        // generate_the_Monomials_matrix(AA,num_sp,spx,spy,spz,c);

        // double *GINVA=new double[num_sp*num_sp+1],

        // double TOL=1e-10;

        // int IRANK;

        generate_the_Monomials_matrix_fortran_format(SA, num_sp, spx, spy, spz, c, _k, lambda);

        /*std::cout << "\nSA:" << "\n";
        int SA_size = num_sp * num_sp + 1;
        for (int i = 0; i < SA_size; i++) {
            std::cout << std::setprecision(16) << SA[i] << "  ";
            if (i % num_sp == 0) std::cout << "\n";
        }*/

        /*   k=0;

           for(i=0;i<num_sp;i++)

           {         for(j=0;j<num_sp;j++)

                    {

                      fot<<"A"<<k<<"="<<SA[i*num_sp+(j+1)]<<"\t";

                      k++;

                    }

                    fot<<endl;

           } */

        for (i = 0; i < num_sp; i++)

            for (j = 0; j < num_sp; j++)
            {

                svdA[(i + 1) * (num_sp + 1) + j + 1] = SA[j * num_sp + (i + 1)];
            }

        // add_c_ToDiagonal(svdA, num_sp + 1, c / lambda);

        // //std::cout << "\n\nsvdA:" << "\n";
        // for (int i = 0; i < num_sp; i++) {
        //     //std::cout << "{";
        //     for (int j = 0; j < num_sp; j++) {
        //         A_new[i * num_sp + j] = svdA[(i + 1) * (num_sp + 1) + j + 1];
        //         //std::cout << std::setprecision(16) << svdA[(i + 1) * (num_sp + 1) + j + 1];
        //         //if (j < num_sp - 1) std::cout << ",  ";
        //     }
        //     //std::cout << "},\n";
        // }

        svdcmp(svdA, num_sp, num_sp, svdW, svdV);

        // //std::cout << "\n\nsvd_U:" << "\n";
        // for (int i = 0; i < num_sp; i++) {
        //     for (int j = 0; j < num_sp; j++) {
        //         svd_U[i * num_sp + j] = svdA[(i + 1) * (num_sp + 1) + j + 1];
        //         //std::cout << std::setprecision(16) << svdA[(i + 1) * (num_sp + 1) + j + 1] << "   ";
        //     }
        //     //std::cout << "\n";
        // }
        // //std::cout << "\n";

        // //std::cout << "\n\nsvd_S:" << "\n";
        // for (int i = 0; i < num_sp; i++) {
        //     svd_S[i * num_sp + i] = svdW[i + 1];
        //     //std::cout << std::setprecision(16) << svdW[i + 1] << "   ";
        // }
        // //std::cout << "\n";

        // //std::cout << "\n\nsvd_V:" << "\n";
        // for (int i = 0; i < num_sp; i++) {
        //     for (int j = 0; j < num_sp; j++) {
        //         svd_V[i * num_sp + j] = svdV[(i + 1) * (num_sp + 1) + j + 1];
        //         //std::cout << std::setprecision(16) << svdV[(i + 1) * (num_sp + 1) + j + 1] << "   ";
        //     }
        //     //std::cout << "\n";
        // }
        // //std::cout << "\n";

        // //ʵ��A����ĸ�ԭ A = U * S * V.t()
        // int m = num_sp, n = num_sp, p = num_sp;
        // //std::cout << "Recover A = U * S * V.t()";
        // double* US = new double[m * n];
        // matrixMultiply(svd_U, svd_S, US, m, m, n);
        // double* svd_Vt = new double[n * n];
        // matrixTranspose(svd_V, svd_Vt, n, n);
        // matrixMultiply(US, svd_Vt, Recovered_A, m, n, n);

        /* std::cout << "\nsvd_Vt:" << "\n";
         for (int i = 0; i < num_sp * num_sp; i++) {
             std::cout << std::setprecision(16) << svd_Vt[i] << "   ";
             if ((i + 1) % num_sp == 0) {
                 std::cout << "\n";
             }
         }
         std::cout << "\n";*/

         /*std::cout << "\n\nRecovered_A : \n";
         for (int i = 0; i < num_sp * num_sp; i++) {
             std::cout << std::setprecision(16) << Recovered_A[i] << "   ";
             if ((i + 1) % num_sp == 0) std::cout << "\n";
         }
         std::cout << "\n";*/

         // double Recover_Avergae_error = 0.0;
         // for (int i = 0; i < num_sp * num_sp; i++) {
         //     Recover_Avergae_error += abs(A_new[i] - Recovered_A[i]) / abs(A_new[i]);
         // }
         // Recover_Avergae_error = Recover_Avergae_error / (num_sp * num_sp);
         // if ((SVD_min_error - Recover_Avergae_error) > 1E-20) {
         //     SVD_min_error = Recover_Avergae_error;
         // }
         // if ((Recover_Avergae_error - SVD_max_error) > 1E-20) {
         //     SVD_max_error = Recover_Avergae_error;
         // }

         ////////////////////////////////////////////////

         // Generate the polynomial coefficients.

         ////////////////////////////////////////////////

         //   SVD(&num_sp,&num_sp, A[1],&num_sp,&TOL,&IRANK,&GINVA[1],&num_sp);

         //   for(i=0;i<num_sp;i++)  //Initialize the value to the polynomials.

         //   for(j=0;j<num_sp;j++)

         //   {

         //   	P[i][j]=GINVA[(i+1)+num_sp*j];

         //   }

        for (i = 0; i < num_sp; i++)

            for (j = 0; j < num_sp; j++)
            {

                P[i][j] = 0;

                for (k = 0; k < num_sp; k++)

                    P[i][j] += svdV[(i + 1) * (num_sp + 1) + (k + 1)] * svdA[(j + 1) * (num_sp + 1) + (k + 1)] /
                    svdW[k + 1];
            }

        /*std::cout << "\n\nP:" << "\n";
        for (int i = 0; i < num_sp; i++) {
            for (int j = 0; j < num_sp; j++) {
                std::cout << std::setprecision(16) << P[i][j] << "   ";
            }
            std::cout << "\n";
        }
        std::cout << "\n";*/
        // std::cout << "Recover_Avergae_error = " << std::setprecision(16) << Recover_Avergae_error << "\n";

        ///////////////////////////////////////////////////////////////////////////////

        // Test the accuracy of the interpolation in some uniformly spaced points.

        ///////////////////////////////////////////////////////////////////////////////

        int num_p = numsp_1D + 1; // numsp_1D*2*2; //num_p: the number of test points in one dimension.

        double dl = length / (num_p - 1);

        double xx, yy, zz; // Used to store the coordinates of the current points.

        Complex G0, Ga; // G0: the true value; Ga: the approximated value.

        dvector rf;

        double R, error, RBF_errorsum, LBF_errorsum;

        double errormin = 1e10, errormax = -1e10;
        vector<double> min_error_XYZ, max_error_XYZ;
        vector<double> RBF_error_list;

        RBF_errorsum = 0.0;
        LBF_errorsum = 0.0;

        TestXYZ_list.clear();
        for (i = 0; i < num_p; i++)
        {
            for (j = 0; j < num_p; j++)
            {
                for (k = 0; k < num_p; k++)
                {
                    TestXYZ_list.push_back(vector<double>{i* dl, j* dl, k* dl});
                }
            }
        }

        // int q = 21;
        // if(CalLBF) std::cout << "Individual testing of the " << q << "-th point" << std::endl;
        // double qth_test_x = TestXYZ_list[q][0];
        // double qth_test_y = TestXYZ_list[q][1];
        // double qth_test_z = TestXYZ_list[q][2];

        // dvector rf_qth = dvector(qth_test_x, qth_test_y, qth_test_z);
        // // std::cout << "i = " << i << "   rf_new = (" << rf_new.vect_x << ", " << rf_new.vect_y << ", " << rf_new.vect_z << ")";
        // // std::cout << "   rs_new = (" << rs.vect_x << ", " << rs.vect_y << ", " << rs.vect_z << ")\n";

        // double R_qth = (rs - rf_qth).abs_vector();
        // Complex Real_G_qth = exp(Complex(0, 1) * _k * R_qth) / R_qth;
        // Complex RBF_G_qth = Complex(0.0, 0.0);
        // scattered_point_RBF_interpolation(RBF_G_qth, qth_test_x, qth_test_y, qth_test_z, GG, P, poly, num_sp, c, spx, spy, spz, _k, lambda);
        // double RBF_error_qth = abs(RBF_G_qth - Real_G_qth) / abs(Real_G_qth);
        // if(RBF_min_error_qthPoint - RBF_error_qth > 1E-10) {
        //     RBF_min_error_qthPoint = RBF_error_qth;
        //     RBF_optimal_c_qth = c;
        // }

        // if(CalLBF) {
        //     Complex LBF_G_qth = Complex(0.0, 0.0);
        //     scattered_point_lagrange_interpolation(LBF_G_qth, qth_test_x, qth_test_y, qth_test_z, GG, numsp_1D, zeroX, zeroY, zeroZ);
        //     LBF_error_qthPoint = abs(Real_G_qth - LBF_G_qth) / abs(Real_G_qth);
        //     CalLBF = false;
        // }

        for (int i = 0; i < TestXYZ_list.size(); i++)
        {
            double Test_x = TestXYZ_list[i][0];
            double Test_y = TestXYZ_list[i][1];
            double Test_z = TestXYZ_list[i][2];

            dvector rf_new = dvector(Test_x, Test_y, Test_z);
            // std::cout << "i = " << i << "   rf_new = (" << rf_new.vect_x << ", " << rf_new.vect_y << ", " << rf_new.vect_z << ")";
            // std::cout << "   rs_new = (" << rs.vect_x << ", " << rs.vect_y << ", " << rs.vect_z << ")\n";

            double R_new = (rs - rf_new).abs_vector();
            Complex Real_G = exp(Complex(0, 1) * _k * R_new) / R_new;
            Complex RBF_G = Complex(0.0, 0.0);
            scattered_point_RBF_interpolation(RBF_G, Test_x, Test_y, Test_z, GG, P, poly, num_sp, c, spx, spy, spz, _k, lambda);

            double RBF_error = abs(Real_G - RBF_G) / abs(Real_G);
            RBF_error_list.push_back(RBF_error);
            RBF_errorsum += RBF_error;

            if (errormin - RBF_error > 1E-10)
            {
                errormin = RBF_error;
                min_error_XYZ = vector<double>{ Test_x, Test_y, Test_z };
            }
            if (RBF_error - errormax > 1E-10)
            {
                errormax = RBF_error;
                max_error_XYZ = vector<double>{ Test_x, Test_y, Test_z };
            }

            if (LBF_error_list.size() < TestXYZ_list.size())
            {
                Complex LBF_G = Complex(0.0, 0.0);
                scattered_point_lagrange_interpolation(LBF_G, Test_x, Test_y, Test_z, GG, numsp_1D, zeroX, zeroY, zeroZ);
                double LBF_error = abs(Real_G - LBF_G) / abs(Real_G);
                LBF_error_list.push_back(LBF_error);
                LBF_errorsum += LBF_error;
                if (i == TestXYZ_list.size() - 1)
                {
                    LBF_Average_error = LBF_errorsum / TestXYZ_list.size();
                }
            }
        }

        double RBF_Average_error = RBF_errorsum / TestXYZ_list.size();

        if (Previous_RBF_Average_Error < 0.0)
            Previous_RBF_Average_Error = RBF_Average_error;
        else if (RBF_Average_error - Previous_RBF_Average_Error > 1E-10)
        {
            Error_increase_count++;
            if (Error_increase_count == 10)
            {
                std::cout << "Error_increase_count == 10" << std::endl;
                break;
            }
        }
        else
        {
            Error_increase_count = 0;
        }
        Previous_RBF_Average_Error = RBF_Average_error;

        // if(RBF_min_Average_error - RBF_Average_error > 1E-10) {
        //     RBF_min_Average_error = RBF_Average_error;
        //     RBF_Min_error_list = RBF_error_list;
        //     RBF_optimal_c = c;
        // }

        if ((0.02 - RBF_Average_error) > 1E-10)
        {
            std::cout << "RBF_min_Average_error < 0.02" << std::endl;
            RBF_min_Average_error = RBF_Average_error;
            RBF_Min_error_list = RBF_error_list;
            RBF_optimal_c = c;
            break;
        }

        // for (i = 0; i < num_p; i++)

        //    for (j = 0; j < num_p; j++)

        //        for (k = 0; k < num_p; k++) {

        //            xx = i * dl;
        //            yy = j * dl;
        //            zz = k * dl;
        //            std::cout << "(xx, yy, zz) = " << "(" << xx << ", " << yy << ", " << zz << ",)\n";

        //            rf = dvector(xx, yy, zz);

        //            R = (rs - rf).abs_vector();

        //            G0 = exp(Complex(0, 1) * _k * R) / R;

        //            scattered_point_RBF_interpolation(Ga, xx, yy, zz, GG, P, poly, num_sp, c, spx, spy, spz, _k,
        //                lambda);

        //            // scattered_point_lagrange_interpolation(Ga,xx,yy,zz,GG,P,poly,num_sp);

        //            error = abs((G0 - Ga) / G0);
        //            std::cout << "error = " << error << "   G0 = " << G0 << "   Ga = " << Ga << "\n";
        //            std::cout << "----------------------------------------------------\n\n";

        //            if (errormin > error)
        //                errormin = error;

        //            if (errormax < error)
        //                errormax = error;

        //            // cout<<"\n error ="<<error;

        //            errorsum += error;
        //        }

        // error = errorsum / (num_p * num_p * num_p);

        // std::cout << "\n Minimum error =" << errormin << "\n";

        // std::cout << "\n Maximum error =" << errormax << "\n";

        // std::cout << "\n Average error =" << error << "\n";

        // std::cout << "----------------------------------------------------------------\n";

        /* fo.open("error_rad_svd.txt", ios::ate);

         if (error) {

             if (test_type == 3 || test_type == 4)

                 fo << num_sp << " d=" << length / lambda << " c=" << c << " souce position="
                 << (0.5 + nc * dlength) * length << " error=" << error << '\n';

             else

                 fo << num_sp << " d=" << length / lambda << " c= " << c << " error= " << error << '\n';
         }

         fo.close();*/

         //   if((test_type==1||test_type==2)&&constant_c)

         //   {

         //   	break;

         //   }
    }

    Length_OptimalC_list.push_back(make_pair(length, RBF_optimal_c));

    // std::cout << "RBF_min_error_qthPoint = " << RBF_min_error_qthPoint << "   LBF_error_qthPoint = " << LBF_error_qthPoint << "\n";
    // std::cout << "RBF_optimal_c_qth = " << RBF_optimal_c_qth << "\n\n";

    // if(RBF_min_Average_error - 0.02 > 1E-10) {
    std::cout << "RBF_min_Average_error = " << std::setprecision(12) << RBF_min_Average_error << "   ";
    std::cout << "LBF_Average_error = " << std::setprecision(12) << LBF_Average_error << "\n";
    std::cout << "RBF_optimal_c = " << std::setprecision(6) << RBF_optimal_c << "\n\n";
    // for(int i = 0; i < RBF_Min_error_list.size(); i++) {
    //     std::cout << "RBF_Min_error_list[" << i << "] = " << RBF_Min_error_list[i] << "   ";
    //     std::cout << "LBF_error_list[" << i << "] = " << LBF_error_list[i] << "\n";
    // }
    // std::cout << "RBF_optimal_G = " << std::setprecision(18) << RBF_optimal_G << "\n";
    // std::cout << "\nsvd_min_error = " << std::setprecision(18) << SVD_min_error << "   svd_max_error = " << std::setprecision(18) << SVD_max_error << "\n";
    std::cout << "=============================================================\n\n";
    // }

    return RBF_min_Average_error;
    // return LBF_Average_error;
}
